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ABSTRACT 


Two numerical methods of determining optimun 
paths in the problem of Bolza are presented, The basic 
theory of the methods is outlined and their essential 
characteristics are demonstrated by several specific 
applications to ship routing. 

Ihe first method is characterized as one of differ- 
ential corrections employing methods developed by Professor 
Faulkner of the U, S, Naval Postgraduate School. The 
second method is termed a method of gradients employing 
calculation routines developed by Professor Bryson of 
Harvard. On the basis of the applications comparisons are 
made in the areas of ease of applicatior and speed in 
producing results. 

I wish to express my appreciation for the guidance 
and encouragement of Professor Faulkner and for the 
programming of Mary Haynes of the Postgraduate School 


Computer Center. 
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INTRODUCTION 

Ihe problem of obtaining maximum or minimum performance 
in a variety of controllable situations is undoubtedly as old 
as time. Usually the degree of success is determined by the 
efficiency with which some vital, limited quantity is ex- 
pended. 

Often the best solution cannot be defined, except as it 
develops by experience, even though the problem ean be de- 
fined and stated mathematically. The magnitude of the 
analysis involved limits the feasibility of determining the 
solution. Modern high speed computing machinery has, for 
many problems, removed the latter limitation. As a result 
renewed interest has been shown and much progress made in 
the development and application of control theory to such 
problems of optimum control. 

Ihis paper presents an outline of two practical nu- 
merical methods of solution of some general problems of this 
type. The class of problems is more specifically termed the 
problem of Bolza [1]. For purposes of identification in this 
paper the two methods will be called the ‘differential 
method! and the method of the fundamental lemma' which will 
henceforth be abbreviated the D-method and the FL-method 
respectively. 

Ihe D-method is based on the formulas for differentials 
due to Bliss [2]. The numerical scheme for determining the con- 


stants of the solution by a Newton iteration is the work of 
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Faulkner [3],[4]. His research papers and related studies 
are the source of both the theory and the numerical data 
presented for the D-method. 

In the literature the FL-method is also referred to as 
the 'gradient' and the 'steepest ascent' method. The calcu- 
lational procedure is similar to the proof of the fundamen- 
tal lemma and is apparently conceptually due to Courant [5]. 
The numerical scheme, which makes possible its application 
to the class of problems being considered, is the work of 
Bryson and Denham [6]. A similar method is presented by 
Kelly [7]. The theory and basic procedure for application of 
the FL-method, as reported in this paper, is that of Bryson 
and Denham. The contribution of the writer lies in the de- 
velopment of specific convergence schemes, left undefined by 
the authors [6], which make possible the programming and so- 
lution of two relatively simple ship routing problems thereby 
demonstrating several features of the method. Both problems 
have also been solved by the D-method and serve as a basis 


for limited comparison of the methods, 
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I 
INTRODUCTION TO PROBLEMS IN OPTIMUM CONTROL 


1.1 Formulation of a problem 
We desire the solution to the set of ordinary nonlinear 


differential equations 


€ 
(1) x, = f,[x.(t),x,(t),....,x (t),p.(t),..,p (t),t]J 


í = 2و۱‎ t. X ۳ 


which is optimal, in the sense of maximizing the terminal 


value of a performance criterion 


ar 
gx... Хь дт > 


while simultaneously satisfying terminal constraints 
(2) һ(х,,......,х , ©), 503 k=1,2,...,M3 O (N (n. 


To effect these conditions we have at our disposal m con- 
trol or driving functions Pg), S- 1,2,....,m. Their 


definition for t. < t < T is a major portion of the problem. 


* Adot ( ) over a variable is used to designate the deriv- 
ative with respect to t. 


xx If an integral is to be maximized we introduce an addi- 
tional state variable, X +]? and an additional differen- 
tial equation 


X. +1 n+ د‎ е е е X t) 
where +] is the integrand of the integral. 
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1 با‎ is assumed the functions f. are defined over an 
open region and are class C", It is also assumed the func- 
tions g and hy. are independent and of class C' in the 
terminal region. The variables (X X5, 9X4) will be 
termed the state variables. All are functions of the in- 
dependent variable t, The terminal value, T, of t is 
either specified or otherwise determined from a terminal 
condition. The values (ху,хь,* бух to are assumed specified. 

It will be convenient to represent these sets of varia- 
bles by matrices or by vectors. We define 


x4 (t) pi CU) ha (x, t) f(x, p, t) 
х = xot) ;p-7|p5,(0)| y h- |h59(x,t)| ; f = fo(x,p,t). 


е 


Xy (4) | Palt) | hy(x,t) | | f£, (x, p, t) 


No distinction will generally be made between an nx] 
(column) or a іхп (row) matrix and a vector. If a quan- 
tity is best considered a vector, as in a scalar product, 
it will be designated, for example, Р. 

Equations (1) апа (2) then become 


(3) fe = “fc pet) 
(4) БОЕ раны О. 


We desire the solution of system (3) which maximizes 
the performance quantity g(x,t) т subject to the con- 
straints (4). For purposes of introduction of terminology 
and some basic theory a brief review of the related elements 
of the Calculus of Variations is first presented, 


1.2 The Classical Calculus of Variations 
The problem formulated in section (1.1) is typical of 
those considered in the variational calculus - those of 
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determining maxima or minima or, in general, staticnary 
values (defined on page 8 ) of functicnals such as g(x,t), 
by seeking the argument functicns x(t} and p(t) for 
which the functional assumes the stationary value or ex- 
tremum in question. 

An outline of scme phases of the classical method of 
the variational calculus, as. applied to a greatly simpli- 
fied optimum problem, will serve to introduce terminology 
and to present essential concepts which are the basis for 
both methods of solution being considered. The method may 
be termed 'indirect' in that it is based on the reduction 
of the variational problem to differential equations the so- 
lution of which, subject to the boundary conditions, is the 
solution to the problem. The two methods considered in 
this study may be termed 'direct' in that they consist of 
construction of a sequence of functions that converges to 
the desired function, 

Following Courent and Hilbert [8], we consider the 
Simplest problem of the variational calculus. Find the 
function y(t) which is such that y(t.) = А апа 
y(t) = B, Fig. 1, and the functional 


T 
FLyCt)] - | f[t,y(t),y'"(t)] at 


0 
is maximized. Any function y = y(t) which meets the end 


concitions, is of class C, and has a piecewise continuous 
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Fiy(t)].  Sugpose we have 


found the desired extremal function y = y(t). Let us con- 
sider siso an unspecified new function n(t) which is de- 


fined for ў $ ts t possesses the necessary derivatives 


O 1۳ 
and continuity, and is such that n(t,) = n(t,) = 0, but is 


otherwise arbitrary. The function 
y = y(t) +t en(t) = y(t) + by(t), 


where e is 2 parameter, represents a one-parameter family 
of admissible functions which contains the solution y(t) 


when e=0, The quantity 6éy(t) = en(t) is the variation 





of the function y = y(t). For sufficiently small magni- 
sudes of e the varied functions y(t) lie in an arbi- 


trorily small neighborhood of the extremal, y = y(t), and 
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the Euler equation. A solution to the Buler egustion 
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called an extremal. Its validity is a necessary 06۳0۵1 
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for the existance of an extremum,[8] page 1855 Being a 
differential equation of the second order, two arbitrary 
constants must be determined for a solution which meets the 
bcundary conditions, i.e., to make an extremal an admissible 
extremal. Every solution of Euler's equation is an ex- 
tremal of the maximum problem,[8]. Courant describes (6) 
the variational derivative of F with respect to t and 
terms its role here as analogous to that of the gradient in 
ordinary maximum problems. 


ln the n+l dimensional case, with coordinates 
ty, Yalt), Yalt), 4... , y Q2 


and the functional 
i | 

(8) I - f(t, Vy» You eves و۲ وی‎ Yo». .» Yp) dt, 
0 

the Euler equations take the form 


۷ 


o E Е a x = 


a system of the second order equal in number to the number 
of unknown functions y, CU to be defined. Апу п+1 space 


C1 NAE 
a nm, Ус” yo(t), ۰ و‎ ۰ ۰ ۰ ө 9 Yn” y, (8) 


which is a solution to the system (9) is an extremal and 


furnishes a stationary value of the functional (8). 
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Ihe necessary condition of Euler is usually called the 
first necessary condition. There are also necessary con- 
ditions of Weierstrass, Legendre and Jacobi as presented, for 
example, in [1] chapter 1. These are not discussed here; 
the fundamental problem confronting an applied mathematician 


is that of finding a likely candidate for a solution, 


1,3 Introduction to direct numerical methods 

The problem as stated in section (1.1) is similar to 
that resulting in system (9). The state variables x(t) of 
our functional g(x,t) are defined by a set of n first-order 
differential equations and one or more control functions, or 
parameters, p(t); t is the independent variable. Solu- 
tions to the related Euler equations will, to the extent of 
meeting the first necessary condition, define the functions 
p(t) such that the resulting path x(t) = 0 “Я ЧІ impart 
at least a stationary value to the functional g(x,t)p. 
Ihe proper choice of the arbitrary constants of these so- 
lutions will result in admissible extremals if they are 
chosen such that the terminal conditions h(x,t)p = O will 
also be met. It will be assumed that the initial conditions 
۰ are specified constants. 
Various indirect approaches to such a solution exist. 
A bibliography to indirect method studies is presented in 
the survey paper [10]. Analytical solutions of the Euler 


equations involved usually require idealizing assumpticns 
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which limit their applicability in practical situations, 
Under more realistic assumptions a numerical attack on 
these equations is required. Some direct methods are pre- 
sented in [8] page 174. In the following sections the 
essentials of two additional practical methods are out- 
lined. As with all methods of solving Bolza's problem, 
Lagrange multipliers are employed producing equations 
equivalent to the Euler equations in terms of solutions 
to an adjoint system of differential equations, The con- 
cept was developed by Bliss [2] in connection with 
differential corrections in ballistics, 

The D-method uses only extremals, Then, by an iter- 
ative correction procedure,obtains an admissible extremal 
out of the family of extremals. 

The FL-method produces an extremal out of any likely 
soluticn to the original equations; part of the problem is 
in a subroutine to approach an admissible solution. Ihe 
solutions to the Euler equations, producing the extremal, 
are approached iteratively by a gradient, or steepest- 
ascent technique, 

Since all methods of solution involve solutions to 
the adjoint equaticns this concept will be outlined be- 
fore the details of the two methods are presented seper- 
ately, 

1.4 The adjoint equations 


The problem as stated in section (1.1) was to determine 
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the solution to the system 


(3) X = f(x y) 
such that the terminal constraints 
(4) h(x,t)4 = O 


are satisfied and the functional g(x, te is maximized. 

Of primary interest is the relationship between vari- 
ations 5p(t) in the control variables p(t), at any 
value t, and the resulting terminal variations ۵ of 
the state variables. We also need the relationships 
between p(t), 6h, and ۰ 

Looking at (3) in component form,(1), the variational 


relationships for the state variables can be expressed as 


AS: 


š af af 
6х. = 15х. + ==1 6х Өрт 


2х1 1 ax Xo 2* "T ах 6x, + m 


i=l, 2, 000 9 Ne 
In matrix notation this is 


(10) 6X = F 6x + G бр 
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where 


е е е = 
8 X] E 
F(t) - : 
f of 
а 8p, 





sp = | X ; бХ = 
6p, Ct) x 6x, 


We will use the prime symbol to denote the transpose of a 
matrix, e.g., 6x! = (6X), 5x5, ....., 6X). 

Following Bliss [2] the system of linear equations 
adjoint to (3) is defined to be 


(12) = -F'(t)A 
where 


AG) 7 D4UO , Ag CD Lee ee (2 


is a matrix of n Lagrange multipliers, and 


í ° ° ° 
NE 1 hy? ho» ۰ 6 ۰ 6 ۰ ۰ و‎ An ilə 


The solutions of systems (3) and (12) are related by 


12 


І. 
‚е 199 


$ =s 


d 13 
EI зе 


Y, gm, <= a” „ J=- -e | 
I uo! ier tae 9799 "ng 


- á ur DS AN ne’ 
& IM LII Jj 


— м ie 

a 
' qi во“ 

La a у Вашей 
ud ona et 





d ; 
dt[(Aex] » AG 6p; 


hence * 
! 1 | ; 
(is) [Asx J, = (AG sp) dt. 
© Ў 
Since we are assuming x, fixed 6x, - O іп (13) and we have 
0 0 
A ۱ T ' 
(14) Их = Fa АС бр) dt 


This is the differential formula relating the terminal changes 
in the state variables and the variations in the control 
variables. It is termed the Fundamental Formula by Bliss [2]. 
It is alsc known as the Fundamental Adjoint Formula and as 


the one-dimensional form of Green's Theorem [11], [12]. 


1.5 Applications of the Fundamental Adjoint Formula 

some indication of the usefulness of the fundamental 
formula is displayed in the following observations. 

1) If A is chosen as the specific solution to the 
adjoint system (12) such that at t = T 


Ai н ш ше о 


RT refers to the gradient of x] in X space, 
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We assume no corners, [9] page 8, exist.‏ 
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let us call this solution ۸ cae (14) then becomes 
Гл 
6x4 CT) - / 6 ۵9 0 و‎ 
O 


yielding the terminal variation of Xy due to variations 
—À— 


6p(t). Similarly choosing A ,(t) such that 


کے 


№, ст) = Уве, (Т) = (030800 ,13+.».2:30) 


i suc И 


gives 


Га 
(15) ёх, (Т) = «Ni Gépdt. 


c) ТЕ N is chosen as the solution to the adjoint 


system such that А.С = Ve(T), (14) implies 
T 


= / 
(16) о) = (УЕ PIE IM G 6p dt 


giving the first variation of the performance function in 


terms of the parameter variations for toS от. 
3) For the constraint functions we likewise have 


І, 
(17) ly P) = IN, نات ۵ تا‎ ў uc ае а 
J t 3 
O 
4) If we were to ignore the problem of meeting con- 


straints and look for a solution which gives a stationary 
value for g(T), we would require the vanishing of the first 


variation,(see appendix of [4]) 


T 
(18) 6g(T) = IN, G 8p dt 
O 
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which requires A 8 =" IAS 


AAT م‎ ei. е е е е е е е af, = O 
ёрї m 
. | 
ef. 0 е е е е е е ef, 
ap} дра 
which represents the m equations 
Ti 
/ Q ۲ абе” 
(19) да * өр = O, 1 > 1...2: ملاو :>> ټك‎ 


i 

Equations (19) with the adjoint differential equations 
(12) constitute what is generally termed the Euler-Lagrenge 
equations [4], [13], 

In equations (18) the matrix AC is clearly an in- 
fluence functicn giving for any time t the effect that a 
unit impulse ép, at that time,wouid have on the perfcru- 
ance function g(T). In that regerd if we write A © as 


& vector, say At), (18) then can be written 


jt 
(тї) = J A - Fb at 
0 


and for maximum change in the performance quantity ЕСТ) 

We would cheese @ parallel to Alt) for СОЕ 

It should be noted that Alt) and &p( t) are m-dimensional 
vectors in psrameter space (p). Courant, [8] page 223, re- 
ters te vector A(t) as the gradient of the functional 


plt) іс Function space. When the optimum soluticn or path 


5 





in function space is obtained this gradient will vanish. 

In the case of our general problem we have, in addition 
to the requirement of maximizing the functioral E(x, ms the 
necessity for meeting terminal constraints (Хх, 1= O. Ме 
will chcose the vector p(t) in parameter space such that 
the path meeting the constraints yields the largest possible 
value for g(T). We now consider two methods of solving the 
problem. 

We have called Е(х,®)т a functional. It is a function 
which becomes a functional since (x,t) represents the endpoint 
of a curve,satisfying the given differential equations with 
jnitial point x( t.) given. The function g has become a 
functional in the sense that the point (x,t) depends upon 


the curve, 
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sl The differential method 


Recalling the fundamental differential formula 
۱ # T / 
(14) [Aéx ]« EA AG брае, 
о 


the D-method is based on the Euler equation which states 
that if the path is to furnish a maximum to some functional 
E(x, Jm, it is necessary that ۱6 189 ۲ ۰ 61 610 با‎ A8 of 6p 

in (14) vanish for some solution A of the adjoint system 
for all бы <a This reguirement, Ло = O, expanded 

in terms of (11), implies that for all t € t. T there 
exists a Am andi x(t) иса that 


D E 2 
Ac) 85 =0 
(t) . Qf s= O 
(20) А “Зра 
i N ۰ af =æ 
A а, 


The problem becomes: determine A (t) and hence the control 
functions p(t) such that (20), (3), (12) and (4) are satisfied. 
Assuming first that T is specified we proceed to out- 

line the method of determining the 9 matrix. On page 14 


a set of n solutions of the adjoint system (12) was 
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developed such that 


کے ۔ 


RACE, TEA. m. 
The linear combination 
(21) A(t) = — * cm + 0900600099 + Ne 


of these n solutions is also a solution to (12). 

As a first approximation to the admissible extremal 
solution desired we estimate a Set of constants с; for 
(21) and determine the associated extremal by employing 
this A(t) in system (20). This will in general not be 
an admissible extremal as it will not meet the terminal 
constraints (4). 

On the basis of the error in meeting the terminal con- 
ditions; corrections de, for the constants in (21) are 
determined. The resulting improved A(t), when employed 
with (3) to obtain a new solution, will, if the method con- 
verges, give a more nearly admissible extremal. The process 
is repeated until some convergence criterion is satisfied. 

If the terminal T is not specified it too will be es- 
timated for the first Solution and a correction dT will be 
produced with the de, at the end of each iteration. 

A more detailed analysis of the basis and the calcu- 


lation procedure is presented in appendix (A). 
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2.2 The method of the fundamental lemma 

This method is concerned directly with the construction 
of the set of m control variables p(t), for t< t< T, 
which produce the solution to (3) that satisfies the terminal 
conditions (4) and maximizes the performance quantity 
gx, t). Essentially the procedure is to produce any reason- 
able trial solution by the choice of a nominal control pro- 
gram b(t). The resulting solution will in general be 
neither admissible nor an extremal. That is, the con- 
straints (4) are not satisfied and the Euler equations 
(19), a necessary condition, are not satisfied. A better 
solution can therefore be constructed. A variation sp(t) 
is then produced such that the Euler equations are more 
nearly satisfied and the solution's deviation from that 
of an admissible extremal is reduced. This variation ap- 
plied to the nominal p(t) produces a better control pro- 
gram Dt) = ptt) + 6p(t). The procedure might be de- 
scribed as iteratively stepping a control vector p(t), 
defined in parameter space for ES ts T, in the direc- 
tion of steepest ascent toward meeting constraints and 
toward meeting the first necessary conditions of an ex- 
tremal. This direction is that of an and VE in 
parameter space, where we use the subscript (p) to denote 
parameter space. As an admissible solution approaches an 
extremal Vie will tend to zero [6]. The direction of 
6p, for each point of a nominal solution, will be determined 
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by the direction cf tbe gradients at that point, The 
magnitude of зр at caca point will be best chosen pro- 
portional to the magnitude of the gradients at that point 
in parameter space, 

11 appendix B the desired AS correction program 
p(t) is derived and defined for t< t < T as 


( 22) -— 


m (43) - ahi Ah -1 
ép(t) =46 Non Ano inh Ing’ I -I/ 171 I а уо Inh ur 
88 hg hh he 


Ihe plus sign is used if g(x,t)4 is being maximized, 

One phase of the calculation procedure as outlined by 
Bryson [6], and summarized in appendix B, was found to be 
scmewhat indefinite, It is that of determining the arbi- 
trary constant (as)? which governs the magnitude of the 
step size p(t) by limiting the integral of these mag- 
nitudes for ES t< T by the relationship (B-7). It 
seems apparent that the larger the magnitude (а=), within 
the limit which allows the linearization (10), the more 
rapidly our solution process will converge to the admissible 
extremal. Accordingly the first restraint on (45) м111 
be imposing an upper limit (as )% as estimated from the 
nature of the problem in conformance with (10). 

Further it is obvious that repeated application of 
large magnitude variations  sp(t) will not allow fine ad- 


ju: tments towzrd the admissible extremal, as the solution 





is approached, unless the magnitude (45) is further mna- 
nipulated, We therefore desire an automatic scheme which 
will rapidly step the nominal solution toward the ultimate 
solution but will adjust itself to finer increments in the 
vicinity of the admissible extremal, The following dis- 
cussion of Bryson's procedure is to develop a basis for the 
design of such an automatic convergence scheme. 

The expression (22) defining the vector зр consists 
cf two components which we define TN and бр. We 


Е 
designate 


(23) bp, = G' Ay. Inn dh 


since it is essentially a vector in p space whose magni- 
tude is proportional to the magnitude of the error in meet- 
ing the constraints resulting from the nominal soluticn. 
Sp, has the direction, essentially, of Ура, Accordingly 
Ру is a component of р which produces a more nearly 
admissible solution. The desired correction toward an 
aámissible solution may be large and may require larger 
values of m than the upper limit, previously imposed, 
will allow. We give construction of an admissible so- 
lution first priority, however, and up tc the limit (ds) 
demand that р be designed to meet the constraints, 

The first term of (22), involving the radical, is 
essentially a vector in p-space whose direction is normal 
to V_h.and, as nearly as meeting constraints will allow, is 


р 
parallel to V8: It specifies a component бр» which causes 


el 





progression toward maximization of g(x,t)p without dis- 
turbing the progress toward an edmissible solution made by 


IN The term 


-l ah 


(45)2- аһ. іта 


is what remains of magnitude (45) 2 for use as magnitude for 


De Z 
a vector normal to EP. If бру demands all of (ds) 


this radical is zero and 5р = Py: Progress will be made 


toward an admissible solution and the functional g(x, t) ep 
will increase or decrease according to whether or not "Sp, 


has a component in or opposite to the direction of V_g. 


p 

Accordingly, the scheme for the choice of the magni- 
tude (ds)^ will be to first provide; within the upper 
limit, a magnitude 


-1 


| “Ара Інь dF | 


for attaining an admissible aan and second, provide 
for a component toward the extremal, adjusted according to 
the estimated deviation from an extremal. We define this 
latter magnitude (A), the magnitude of E We then re- 


quire 


(24) (as)? > (as) گم‎ + [6* A, y 17, dh | 


The initial choice and automatic adjustment of the 
magnitude (A) is another problem. We want it as large as 


possible initially. This will be controlled by (24). 
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se want smaller magnitudes as the ultimate scluticn js 
appreachec, automatically adjusted to the need. The con- 
vergence criterion suggested by Bryson [6] will now be em- 
ployed. He develops the fact that on an admissible solu- 
tion the denominator of the radical of (22) is the direc- 


tional derivative or gradient 


(25) s =| Tes تسین‎ Le 
in p-space. This must decrease as the extremal is approached. 
This quantity is, however, in general greater than zero. It 
is reasoned that if a relatively large (A) is repeatedly 
applied to an admissible solution it will rapidly decrease 
the quantity dg/ds of (25) and then will ultimately cause it tc 
pass through a minimum, greater than or equal to zero. If 
at this point the magnitude (A) is not reduced subsequent 
solutions will be over-corrected toward an extremal; trajec- 
tories which oscillate about the solution to the problem. 
we use this characteristic to trigger a reduction of (A) 
to some fraction, say O.2 A, and step from there toward 
the extremal in smaller steps. The minimum value for dg/ds 
attained this time will be smaller and the same behavior will 
signal the need for further reduction in (A). The precess is 
repeated until some lower limit on (A) or dg/ds, a conver- 
gence criterion, is attained, 

A flow chart of a convergence scheme is presented as 


appendix C. 
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в} 
APPLICATIONS AND COMPARISON 


3.1 A basis for comparison 
The performance cf both methods was investigated for 
two ship routing problems, Of primary interest were: 

(a) the problem of choice of starting para- 
meters that result in convergence to a 
solution. 

(b) the domain of these parameters that pro- 
duces a solution. 

(c) the range of problem variaticns that can 
be solved with an acceptable set of these 
parameters, 

(d) the rate cf convergence to a solution. 

For each application the performance criterion was 
defined as the functional 


(x - t/100) + Ay* 
T Е 2 
(26) g = [ с+е ] 6 
ل‎ 


which could be interpreted, for example, as a probability 
of detection cr as a measure of accumulated radicactive 
fallout in an area where the intensity is defined by the 
time and position function g(x,t). This is essentially 
the problem considered by Faulkner [4] and Cook [14]. 


Their programs fcr the D-method were employed as a source 
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of comparison data for that method. 

More specifically the problems considered are: 

Problem I: Determine the route from a point within an 
area where g(x,t) is defined, to any point 
at which g(x,t) is a prescribed tolerable 
maxjmum, such that the accumulated fallout 
g(x,t), is minimized. | 

Problem II: Determine the route between two prescribed 
points in the area where g(x,t) is defined 
such that g(x,t), is minimized. 

For the D-method the calculation procedure of page 18 
for these problems is reduced to the specification of a 
pair of constants defined by Cook [14] as (Ом) and are 
the two components of the initial adjoint vector, which are 
unknown at the outset of computation. 

For the FL-method the starting parameters are, as for 
any problem, a nominal p*(t) and a magnitude (85). For 
these problems p*(t) is a bearing angle or heading and 
could be any constant direction p*(t) =k radians, rough- 
ly away from the point of maximum e(x,t), "ground zero'. 
From (B-7) an acceptable value for (45) is determined by 
(às )* = (8p) “T, where T is approximately the total time 
required and 6p is the change in bearing, at any point, 
that would reasonably meet the linearity requirements of 
(10). As an example; estimating T as 100 with an allow- 


able 5p of 0,5 radians would give: 
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(95) = (0.5) (100) = (5) 


3.2 Problem I 

we consider the specific case of problem I with start- 
ing point (х,у) specified as (0.8,C.1). For the FL- 
method the choice as = 5 is reascnable and a logical 
choice for p*(t) is p*t) = 1/2 for 0S ts T or very 
roughly 'away' from g(max) at (0,0). For the D-method a 
reascnable choice (with the advantage of hind-sight) for 
(ХЧ) is (0,-100) which is a vector reughly in the direc- 
{ісп оѓ the negative gradient of g(x,t) with a magnitude 
approximating the final time Т. 

In an effort to experimentally determine the intervals 
of convergence an interval of values containing these esti- 
mated starting parameters was explored for both methods, 

The results are summarized in Table I and Figures 2 ana 3. 
For the FL-metnod the process converged Гог 

0.1 < равоў 3,2 radians, The starting parameter is 
therefcre any bearing angle rcughly 'away' from (0,0) within 
about a 180° interval. In Fig. 2 path A depicts the op- 
timum route. For values in the vicinity of 31/2 the local 
minimum at the opposite side of the cioud (path B in the 
sketch) is developed. If headed in the direction pact) = 1۲ و‎ 
however, the local maximum will not be produced but the 
nearest minimum (path A) will be generated, while p*(t)= 3,6 


results in convergence to path B, 
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Table i Starting parameters for problem 1. 


EE tee a | "DT ree 





D-Method FL-Method 
NON HO Convergence p (t) Ga er CC 
90 -10 yes O. 1 yes 
80 -20 7 yes 0.2 yes 
60 -40 yes 0.4 yes 
40 -60 yes 0.8 yes 
20 -80 yes a yes 
0 -100 yes раю yes 
= 20 -80 no c O yes 
-40 -60 no 2.4 yes 
-60 -40 no 206 yes 
-80 -20 по 3.0 yes 
-90 -10 no s yes 
0 -50 no 226 (path B) 
0 -200 yes 
60 40 (path B) 
80 20 " 
0 100 ü 
-60 40 no 
-80 20 no 
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converge to path B 


Figure 2 FL Convergence Interval for p (t)=k 


Ihe number of iterations required for a 'solution' de- 
pends on the accuracy demanded. For values within the inter- 
val of convergence the poorest choice for p*(t) requires 
about ten percent more iterations than required of the 
best constant choice. No great penalty is associated with 
rough estimation of the starting parameter,for a close 
approximation to an admissible route is produced rapidly. 
Computer time is of the order of 45 seconds for 40 itera- 
tions, The first 15 iterations produced a very accurate 
soluticn. 

For the D-method starting parameters were explored in 


the intervals -100< 2, < 100, -100 < ug < lOO. For pur- 


0 
poses of probing the influence of magnitude ) و200 (ی لاو ی‎ -50( 
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End (C,-9CC) were «lso tried, Figure 3 sno Table 1 sun- 


мё 


BEIDE cae PERA Ts, іп each case ТОИ СПЕ Оі са 
existed convergence was relatively ravid requiring about 


lO iterations and computer time in the order of 10 seconds 


no Solution t th B 
сё 6,100) A converge to path 
i 
y 


Й 








N 
m COED 


x 
(SO, -10) 


| converge to 
path A 


Figure 3 D-method Convergence Interval for (A suo) 


for the particular convergence criterion employed. 
Continuing with problem I we explore the domain of 
starting points (x Yo) for which a given set of starting 
parameters produces a solution, For practical purposes it 
would be desirable that a given set of parameters func- 
tion for a large domain of starting points under the radio- 
active cloud, Accordingly solutions were attempted for a 
number of starting points in an interval bounded roughly by 


the left and right extremeties of the intolerable fallout 


eo 






А baa ( ПЕЕ o Е 


о ў а» лы - д! 
ЕЕ > oS 1 7 


area. The input parameters were maintained at (15549) 7(0, -100) 
for the D-method, and p*(t)=1.0 for the FL-method. 

Figure 4 and table II portray the results for eight start- 
ing points, Both methods produced solutions for each run, In 
figure 4 the terminal points are associated with the corre- 
sponding initial point by a dashed line simulating the route, 
As table II indicates, the coordinates of the terminal points 
agree for the two methods to the third or fourth significant 
figure. Values for final time (T) and the performance g like- 
wise agreed to at least four figures, Extended results for them 
are given only for run Рр, 

In regards to the FL-method, table III is a printout of 
the terminal values of the iterations toward a solution for 
run D of figure 4 and table II. It will serve to display 
several features of the FL-method as applied to a relatively 
extreme choice of the starting parameter p*(t). Specifically 
the starting point is (0.3,0.1) and p*(t) was, for probing 
purposes, chosen as p*(t)=0.1. This is an obviously unwise. 
choice in as much as the cloud and the ship have comparable 
Speeds and, for this case, are almost parallel. Consequently 
the terminal time (T) for the nominal run is about ten 
times the value on the extremal, A nominal solution is pro- 
duced, however, and from it relatively rapid convergence to 
an extremal is portrayed. 


The parameter (ds) was chosen as (6.0). The choice of 
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Fig. 4 Allowable Starting Points for (1,,u,) and p (t) fixed 


Table II Allowable Starting Points (x Yo) for Problem 1 
x 
(> ue) - (0851000; TES = 1.0 


Start point X(T) | X 
(x Yo) 1 | A : 

À (1.2,.1) | 1.4092 1.4103 | ICE 170516 
Boe. om) x 1.1379 1.1388 | 1.0702 1.0702 
С (в) | 0.7266 | 0.7223 | 1.0661 1.0662 
D (ommo | OS x 0.1778 | 1.0074 1 б С 
E (-0.7,0.1) -1.1066 x -1.1071 | 0.6128 0.6125 
Е 1) «174175 | -1.4179 | 0.4406 | 0.4401 
С --0.1) | 1.5879 | -1.5882 0.2325 | 0.3320 
H (-1.4,0.1) | -1.7292 x -1.7293 0.2454 | 0.2452 


Typical values for final time T and performance g(T) 


for point D: T,, = 91.6399, Tp = 91.6417 
в(Т) = 49.0500, G(T)p = 49.0501 
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Xo = 0.3000 ү 


d 


967.1896 
124.8307 
95,2305 
98.7530 
94,0974 


91.8374 
92.2926 
91.8649 
91.6341 
91.7693 
91.6957 


91.6437 
91.6199 
91.6262 
91.6340 
91.6426 
91.6293 


91.6331 
91.6372 
91.6416 
91.6376 
91.6388 
91.6399 


91.6395 
91.6396 
91.6397 
91.6399 
91.6398 
91.6399 


91.6399 
91.6398 
91.6399 
91,6399 
91.6399 
91.6399 


Computer time: 


Table III Problem I, FL-Method, Run D 


0 
X(T) 


9.9236 
2.6762 
1.0787 
.3507 
. 1648 
„1323 


.1091 
© 2244 
. 2062 
. 1808 
„1826 
„1801 


AO 
.1759 
„1765 
.1771 
„1776 
. 1766 


«1769 
“1773 
„1776 
779 
1774 
1775 


.1774 
.1774 
«1774 
1774 
1770 
.1775 


“1775 
“1775 
“1775 
«1775 
„1775 
TO 


= 0.1000 


Y (lo) 


1.0656 
1.0716 
1.0696 
1.0300 
.9910 
. 9939 


. 9938 
1.0146 
1.0727 
1.0081 
1.0078 


1.0074 
1.0072 
1.0073 
1.0074 
1.0074 
1.0073 


1.0073 
1.0074 
1.0074 
1.0074 
1.0074 
1.0074 


1.0074 
1,0074 
1.0074 
1.0074 
1.0074 
1.0074 


1.0074 
1.0074 
140074 
1.0074 
1.0074 
1.0074 


44 seconds 


S = 6,00 


g(X,t)4 


487.8577 
127.6031 


62.3428 
49,8881 
53.9286 
50.3492 


49,1760 
49,341] 
49.1077 
49,0658 
49.1008 
49.0645 


49 53 
49.0569 
49 ¿0528 
49.0505 
49,0501 
49.0503 


49.0501 
49.0500 
49.0501 
49.0500 
49 ¿0500 
49 ¿05900 


49,0500 
49.0500 
49.0500 
49 0 
49.0500 
49.0500 


49.0500 
49.0500 
49.0500 
49.0500 
49.0500 
49.0500 
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p (t) = 0.10 

D ds 
17661.3315 6,0000 
2228.5757 6.0000 
487.7631 6.0000 
79.2727 6.0000 
275.8876 1.8000 
132.3472 1.8000 
29,3284 1.8000 
64,2659 o 5400 
23.8514 . 9400 
15.8866 . 5400 
28.9107 „1620 
15.1434 „1620 
1,5992 «1620 
11.0875 ‚0486 
7,1449 ¿0486 
3.2402 ‚0486 
. 7613 ‚0486 
3.1283 ‚0146 
1.9116 ‚0146 
. 7058 .0146 
.4999 ‚0146 
. 7230 .0044 
. 3542 ‚0044 
.0208 .0044 
‚3385 ¿COLS 
‚2254 ‚0013 
@ш су ‚0013 
. 0149 ‚0013 
‚0327 ‚0004 
‚0116 . 0004 
‚0062 «0004 
‚0189 ‚0001 
‚0106 «0001 
"023 .0001 
‚0040 ‚0000 
‚0000 





any number O < ds < 20 would also produce solutions. Actually 
values nearer 20 would hasten the convergence while values 
smaller than six would require a greater number of iterations. 
Table III also portrays the operation of the convergence 
scheme discussed on page 23. Since in this case there are no 
constraints, A in (24) is (ds). As table III indicates, dS 
stays at its initially chosen value of (6.0) for four itera- 
tions. By this time the extremal has been closely approxi- 
mated not withstanding the difficult nominal route chosen. For 
the fifth iteration dg/ds has crossed through its minimum and 
the reduction of (ds) to 0.3(dS) has been switched. The pro- 
cess then continued until an extremal was crossed again,caus- 
ing a further reduction in step size. In this case the pro- 


4 


cess was halted when the condition dg/ds < 10 ‘was attained. 


3.3 Problem II 

The two methods were less directly compared for problem II. 
Experimental probing was directed primarily at weaknesses of 
the convergence scheme devised for the FL-method. The first 
specific problem of type II was: Determine the optimum route 
from (x5,y9)7(0.03,0.1) to (Xe53 9) 7(1.0,1,.0) with c 3© 
in (26). An admissible solution, simulated by route (A) in 
figure 5, was closely approximated in the third iteration. 
Convergence to an admissible extremal, however, was not re- 
alized. Longer and longer admissible routes, (B), were devel- 


oped, each better than the previous from the standpoint of 
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decreasing the performance criterion, g . It then became appar- 
ent that the solution to this problem does not exist. The quan- 
tity g clearly would be minimized by escaping from under the 
radioactive cloud to some very distant point, returning after 
the cloud has moved away. The FL-method developed a sequence 
of solutions which approached this solution. The D-method for 
this case returned an'error' printout in as much as one of the 
input parameters is a relatively close estimateof the final 
time T. This value being undifined the process stopped due 
to too large an error in the estimate of the magnitude of this 
parameter, 

Ihe FL-method was then applied to the same problem with 
c 0.1 in (26). An admissible extremal resulted in a reason- 
able number of iterations. 

Ihe importance of the careful choice of the stopping con- 
dition Q(x, t) p= O, for the FL-method, was demonstrated by 
problem II for the route (0.03,0.1) to (0.0,1.0),.figure 6, 
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Figure 5 Problem IIa Figure 6 Problem 6 
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with N(x,t)p= O chosen as (x - Xp)p = O, The resulting 
situation is that the stopping constraint is nearly parallel 

to the route making T, and thus x(T) and ECX, t) particularly 
sensitive to variations in the control variable, Consequently 
difficulty arose in that the constraint (x - Xp) q = O was 
difficult to meet with the accuracy desired, The choice of the 
line (y - yg), = 0, or the circle through (Xg,Yp) with 

center (XL Yo? > would alleviate the situation. In general 

)() و26‎ t) m =0 and gx, t). = О should be as nearly parallel 

as possible. 

The difficulty arising from too large a value for A, was 
demonstrated for the FL-method by the solution of problem II 
for the route (0.03,0.1) to (-1,0.1), simulated by curve (D) 
of figure 6. The initial choices of magnitudes ten and five 
Bor Do prođuced very slow convergence toward an admissible 
route, A, = О,1, however, converged reasonably well. This 
can apparently be attributed to errors due to linearization 
on a path with relatively large curvature, 

The D-method produced rapid convergence to solutions for 


each of the last two cases. 


3.4 Conclusions 

On the basis of this study the following observations 
regarding the points of comparison (a) through (d) of 
section 3.1 seem justified. 
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(a) The choice of acceptable starting parameters is di- 
stinctly easier for the FL-method. The matrix p(t) is the 
set of actual control parameters for the problem aná as such 
it is usually a set of well understood physical quantities. 
These are less abstract than the adjoint vector or the con- 
stants Ci of the D-method. General knowledge about the 
physical concepts involved usually permits estimation of a 
starting control program which relatively closely approximates 
that of the optimum solution. 

(b) The interval of convergence, i.e., the domain of 
starting parameters from which a solution can be produced, is 
Significantly larger for the FL-method. Any reasonable esti- 
mate for the control parameter will start the solution. Table 
I and figure 3 indicate there are seemingly good estimates 
for و (ملاوی()‎ such as (-20,-80) which do not produce solu- 
tions for the D-method. 

(c) Once a workable set of starting parameters is found, 
however, both methods produce solutions from it for a wide 
range of variations of the problem. 

(d) If the D-method eonverges it consistently results in 
the more rapid convergence. For similar convergence criteria 
the computation time ratio is about one to two, and the number 
of iterations required, a ratio of about one to four. The 
criteria employed were different for the two methods, however, 


and it appears from table III that the requirement dg/ds < 1074 
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was too tight for practical purposes, The precision is avail- 
able, however, 

Ihe FL-method took only about ten percent more iterations 
to converge from a relatively wild estimate for p (t) than it 


took for a very close estimate, 


Ihe following comments regarding the two methods are 
based on experience in implementing the two methods and on 
analyses of the papers by Faulkner [4] and Bryson and Denham 
[6]. 

The requirement of storing an entire control program and 
a complete set of functions Агаё for any one iteration is a 
disadvantage for the FL-method . This could result in memory 
storage problems for large problems or small computers, In 
contrast the D-method holds only the parameters С; and T or 
their equivalent. 

As indicated in Appendix B the:FL-method has been develop- 
ed to the point of providing for variations at = and for 
the insertion of weighting functions by the matrix W(t). The 
power of these features has not been investigated in this 
paper. 

The problem of discontinuous control (corners) or un- 
bounded control have also been ignored. The D-method has been 
developed to handle them, [4], while the application of the 


FL-method to such problems has not been undertaken, nor is it 
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clear how this could be done, 

With the D-methoä there is noe question of how closely we 
have approximated an extremal,with the mathematical routine, 
since cnly extremals are used, The only limitation is the 
accuracy of tne integration method. Also, since only extremals 
are used various identities associated with an extremal can 
be employed. For example, if the independent variable does 
not enter explicitly, the Hamiltonian is constant. 

The integral relations for the corrections in the FL- 
method seem to be more stable from the computational standpoint, 
while, if the strong Legendre condition is not satisfied the 
routine must be modified in the D-method. 

Only forward integration (or only backward integration) 
is used in the D-method, In the FL-method the basic equations 
are integrated forward and the corrections are obtained by 
backward integration of the adjoint system. This necessitates 
either storing the complete trajectory just produced or re- 
constructing it from terminal values by backwards integration 


of the basic differential equations of the system. 
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APPENDIX (A) THE DIFFERENTIAL METHOD 


å.l Outline of the analysis 
The following relations have been defined: 


(3) x = f(x,p,t) the state variable equations. 
(4) h(x,t)4 N terminal constraints. 
(12) А а. “ВД the adjoint system, 
T 
(14) [A' 6x]. = " A' G əp dt the fundamental 
to differential formula. 
(21) (0) = сулу + е +... + е, = hi 
e 
^n 


where КО = [р is the solution to (12) such that 
23 


Ani 
NO) 9 31,173 
0,173 
n 
and A. 7 2 саў. 


ј=1 
We write (14), using (21), as 


T 
(A-1) [ ۱۱ 6x], 3 1'(t) G 6p dt 
0 
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From (20) and (21) we get the m equations 
(A-2) A'(t) G(t) = 


It will be recalled that (12) with (A-2) constitute the 
Euler-Lagrange equations. 

Additional relationships between the quantities dT, 
бхр, бр and dc, will now be developed, See [4]. 


A differential dT in T results in the differential 
dx (T) = 6x, (T) + x,dT , i - 1,2,...,n , 


or 


Чхт = бхр + х ат 
giving the variations 
(A-3) бх = (dx - хат). . 
(A-3) with (A-1) gives 

T 

(A-4) \'(dx - x dT)p= [Xàx-Ax ат] =f [X(t)GCt)8p(t) Jat, 
If x is held constant in (A-2), which in expanded form is 

1'G = (e, A+ cohot o EE c A) GC up, t) = 0% 
we can write 


۱ / 
aG " 
m" A.G ас.) +..... AE 6 р + Лб адс) = 


(с A, 


= 
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where we define the m xn matrix 


дб eG eG 
6 == 6 + о е о е е + - ° 


In terms of (21) this is the set of m equations 


OG / ۱‏ ا 
(A-5) À —p + Лус de, + Лб de‏ 


1 
др tert A, б ас, = 0. 


2 


Solving (A-5),for 6p4, 6po,......, бру and substituting 
into (A-4) we get n equations, in the n+l variables 


of the form‏ و وه و۰۰۰۰ وت و1 
n г. Я‏ / 
(А-6) [ ^,ах ZG: + Дух GT Jp , 1 = 1,...09n‏ 


The differential of the constraint n, is 


n 
h = 
> eH dx, Е кат | T 9 k a وه و‎ ot ee e 


( A-7) dh, = ۲ 9x; 
= 


For admissible curves hy is specified, The end point must 


thereforesatisfy the relation 

( A-8) dh, = 0, k = 1,2,....,N. 

In vector notation (A-7) and (A-8) may be expressed by 
(А-9) ( vh, • аг ]. = 03 к= 1,2,..... 30; 


where va, is a gradient in x,t space. 


For a maximum value of E(x, we similarly have 
(A-10) ] ۷۵۰۵۲ Jp = 0 


in x,t space indicating the differential of the terminal 
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point must be such as to result in displacements (geomet- 


tically speaking) normal to the Vg in x,t space, that is, 


tangent to the surface E(x, = 
g(x,t) 


g(maximum), 


Since the functions 


and h(x,t) are assumed 


to be independent the N+1 by n+l matrix of coefficients 


from (A-9) and (A-10) 


oh; eh; о е о е Qh, ah- 

0х1 0X5 Ox, at 
(4-11) ub. 

дру By € oO O O дру ah, 

0x4 Xo 0x, at 

92 o£ обо EE ag 

8X4 9х2 Ox, at 


has rank N + 1. 


The transversal condition, [4] and [2] page 196, may 


be stated as the condition that the rank of the N+e by 
n+1 matrix 
M 
(A-12) M = 
сооро ә ۱ 


be №1 and the rank of the matrices obtained by omitting 


either of the last two rows also be (+1, Row number N+2 


of M is the row of coefficients of dx and dT of (A-4). 
These conditions on (A-12) yield n-N relations involving 

the terminal values of x, 1, and t which must be satis- 
fied by a solution which gives at least a stationary value 


to We designate them 


g( x, t). 
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(A-13) НХ, От = و0‎ S = ое. 


The differential relations derived from (A-13), at the 


terminal point, are 


(кам) ® jam, = 7-04 dx 2.2 Ed ад pa وا‎ + ше) aT 
s at T 


1 9х, i j 0X,0C, aX, ot 
S = ly230...oyn-Ñ, 
where ید‎ is defined as in (21). 
Equations (A-6), (A-9) and (A-14) are 2n equations 
in the “+l differentials dx, ; aT and dC; i and j= 1,.n. 


А.г Calculation procedure 

(a) A solution of the problem is characterized by the n+l 
parameters »یناه ۰۰۰ وچ و دنل و1‎ We arbitrarily choose one | 
relation among the Ci. We might choose AC = 1 or just 
designate one, say Cr” 1. Ме must be careful only that the 
designated Ci does not turn out to be negative or zero, We 
then guess a set of values for the remaining parameters and 
compute simultaneously a solution to the original differential 
equations (3), the adjoint system (12) and the Euler 
equations (A-2) or (20). This yields an extremal Е, 
which does not, in general, satisfy (4) or (A-13). 

(b) Simultaneously calculate the correction integrals (А-6). 
(c) Treat the differential relations (A-5) through (A-13) 


as differences evaluated at the end point of Е. We then have 
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(A-15) An, = h; = h, (t,x)p 5 ай = У. 


AR, 


4 -H, (2), 1 = NEL Ta Ta гае. е Га 


(a) Solve for the n differences و درا شوه ۰۰۰ وتا‎ AT and 
apply to the previous estimate for these parameters. If the 
original estimate was close enough the new parameters will 
result in a more nearly admissible extremal. 

(e) The process is repeated until the sequence of extremals 
converges to an admissible extremal within the limit cf some 
convergence criterion. Additional analysis is required to 


distinguish between maxima, minima and stationary values, 
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Appendix B THE FL-METHOD 


B.1 Outline of analysis 
We have previously defined the relationships 


(3) XE ربا و0 و12‎ n state variable equations 
(4) h(x3t)p = 0 k constraints 
(10) 6x = F 6x G бр variational relationships 
(12) A= "A adjoint differential equations 
T 
4 
(14) ex), / Ле 6p dt fundamental differential 
O formula 


Following Bryson[6] we define one of the functions of 


(4) as a stopping condition, 
(B-1) Ax, t)4 =G 


which produces a terminal value T for t, thus eliminat- 
ing the need for guessing T. 
In section (1.5) we developed, from solutions of the 


adjoint system, 
T 

(1⁄2) sh;( T) = Дм, в 6p dt, „= 1,2; 0 
O J 

or equivalently 


т 
(В-2) &h(T) = TA G 5p dt 
0 
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where A, Ct) is defined as 


Du À е с е © А 
пуу Эһ, Ву 
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| 1 kn eae 


We also developed 
T 
/ 
(8) bet) = fA сорав. 
о 5 


We now add, for (B-1), a similar relationship 


Т 
(B-3) 6n(T) = dA, с бр dt 
о 


where Ao is the solution to (12) such that 


— 


А ст = RN) ° 


where V. denotes a gradient in x-space. 


ror small perterbations the value of T will be changed 


oniy a small amount dT, so that 


dg = ôg + gaT 
Еа) dh = h + ha? 


dn = за + рат 


«here с, п, О, ere defined аз 
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(B-5) o ah . ahr 
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Accordingly (18), (B-2) and (B-3) modified for a change 
dT in T become 


T 
dg = f Apg G 8p dt « ат 
O 


(B-6) T | 
dn = /, A, G sp at + В ат 
O 


T 
an = / Ла 9 өр ае + А ат ۱ 


In order to limit the magnitude of the step ёр 
Such that the linearization (10) is reasonable Bryson im- 


poses the limitation 
2 P ! 
(B-7) (ds) = E op" ¿PIT 


where (dS) is a magnitude arbitrarily chosen as a reasonable 
upper limit to the integral (B-7). 

Throughout this paper a simplified version of the FL- 
method, as developed by Bryson, has been employed in that 
variations at t = b have not been considered and in that 
his matrix of 'weighting functions', W(t), has been 


neglected. 
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sii"incticn ci © jr equations (B-6) gives 


(7 
i 


(b-&) 9 > I A en O OE 
m i 
(B-2) dh = АЛ G &p dt 
where 
е ї E 
М Pe = N - ج‎ ۸۰ ۳۳ , 
gn 5 A п A. h N N 


From (B-7), (8-8) end (B-9) we form the linear 


combination 


T 
(В-10) ава (Aga 9 -v'Ang G - usp) êp dt + v'dh +y (a5)? 
н ۱ 


The second variation Or (B-10) ie 


T 
5006) = ff (AL, G -v'At, G - 2u5p') &^p at 
to en nn EP p 


e 


For maximum dg the coefficient of 6 ^p = © implies 
Bm СА -Л № 
. po nme 


zliminating у and и between (B-7), (B-9) and (B-11) we 


have 


(B-1Z) 





= 


рас (Й. Д. Icl I.) |- ——— ——ÀÓáàhà — E + GA. 17+ an 
ES une go hQ hh hg |, agi glr ho “hh 


(ds)?^- dh'I 


* Use the plus sign for maxima. 
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Calculation procedure for FL-method 

1. State a reasonable nominal control variable program, 
p(t) which is defined for the time period (t 9T). In 
many problems this can be a constant program, For ex- 
ample if the parameters are two directions, say the 
vertical angle | and the bearing ©, we could esti- 
mate p (t)'= ($5,062, where фо is any constant value 
such that O < 0, < 1/2 and ©, is somewhere in the in- 
terval (correct heading + 1/2). Usually better estimates 
can be made, 

2. With p (t) solve system (3) determining T and x(t) 
for to < t < Т. 

3, Solve the adjoint system (12) backwards simulta- 


I... Store 


neously building Aho?’ Ago? Ihh?’ Ing’ gg 


№ y and AG. 

4, On the basis of x(T) from step two, evaluate 
h(x,t)p . Set dh = -h(x,t)m . 

5. Select a reasonable value for dS. If the numerator 
of the radical in (B-12) is negative scale down dh 


such that the numerator is > 0. 
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6. Produce 5p (t) from (B-12) for t. < L لوا‎ 
simultaneously forming p(t) = p (t) + 5p (t) and 
store for the next iteration. 

7, Return to step two forming a new solution until the 


halt criteria are satisfied. 








Appendix C Flow diagram for adjustment of (ds) 
Definition of symbols 
dh -the error in meeting the constraints h(x,t)p = O. 
6 -the magnitude of dh within which it becomes 
reasonable to be concerned about extremality. 


-the convergence critericn for extremality. 


E 
e 
A; -a part of (ds) whose magnitude is adjusted down as 
the solution approaches an extremal, 
D; | -the denominator of the radical of (B-12), 


dg/ds in p-space. 
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